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We present a new numerical method for studying the dynamics of quan- 
tum fluids composed of a Bose-Einstein condensate and a cloud of bosonic or 
fermionic atoms in a mean-field approximation. It combines an explicit time- 
marching algorithm, previously developed for Bose-Einstein condensates in a 
harmonic or optical-lattice potential, with a particle-in-cell MonteCarlo ap- 
proach to the equation of motion for the one-body Wigner distribution func- 
tion in the cold-atom cloud. The method is tested against known analytical 
results on the free expansion of a fermion cloud from a cylindrical harmonic 
trap and is validated by examining how the expansion of the fermionic cloud 
is affected by the simultaneous expansion of a condensate. We then present 
wholly original calculations on a condensate and a thermal cloud inside a har- 
monic well and a superposed optical lattice, by addressing the free expansion 
of the two components and their oscillations under an applied harmonic force. 
These results are discussed in the light of relevant theories and experiments. 
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I. INTRODUCTION 

The physics of ultracold atomic vapours under magnetic or optical confinement has been 
a continuing and ever expanding focus of interest since the realization of Bose-Einstein 
condensation [Q]. Following characterizations of the basic thermodynamic and dynamical 
properties of condensates 0, a number of experiments have been addressed to investigations 
of their phase coherence and superfluidity [H, to the study of non-linear effects and special 
spectroscopic features and to the observation of vortices p. Parallel efforts are being 
made in the study of gases of fermionic atoms and of boson-fermion mixtures [0], with 
the ultimate aim of realizing novel superfluids. 

Theoretical studies of the dynamics of these systems, involving the analytical solution of 
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approximate models have very often been successful in explaining or predicting such novel 
phenomena. However, the interplay of different species or the thermal fluctuations of a 
condensate are not easily handled by analytical methods. The validity of the approximations 
may be limited to the extreme coUisionless or hydrodynamic regimes, the confinement of the 
sample being treated within a local-density approach. Thus, while the equations governing 
these dilute systems remain simple, their numerical solution can be helpful for investigating 
more complex dynamical problems where an intermediate regime is met or a cold-atom cloud 
accompanies the condensate. 

The atomic interactions in a highly dilute Bose gas at very low temperatures, as is 
relevant for current experiments, are described by a contact pseudopotential accounting for 
s-wave scattering in binary atom-atom elastic collisions. A condensate is then treated within 
a mean-field approach by solving the stationary or the time-dependent Gross-Pitaevskii 
equation (GPE). Several types of numerical approaches have been developed: eigenvalue 
solvers [§], variational solvers or explicit solvers fl^ for the ground state, and implicit 
| TT| or explicit |]T2[ time- marching schemes for the dynamics. 

Methods for studying the dynamics of an isolated cloud of ultracold (bosonic or fermionic) 
atoms are also well developed. One needs to solve the Vlasov-Landau equation of motion 
(VLE) for the Wigner distribution function. Various numerical techniques have been used 
for this purpose, which are based either on an ergodic assumption or on the inclusion of 



statistical noise 0], or else use a direct simulation MonteCarlo (DSMC) approach borrowed 
from the methods of molecular dynamics [|T^-[T^. 

The general theoretical background is provided by the book of Kadanoff and Baym 
||T8| . These authors developed a Green's function approach to transport phenomena, which 
extends the Boltzmann equation to strongly interacting quantum fluids and allows for pro- 
gressively improved self-consistent approximations. This formalism was extended to a ho- 
mogeneous Bose-condensed gas at finite temperature by Kane and Kadanoff |]19|, within a 
Beliaev approximation including interactions up to second order. More recently, these meth- 
ods have been adapted to the theory of transport phenomena in a confined Bose gas within 
the Hartree-Fock-Bogoliubov approximation ||20|-p2|, dealing with a Bose- Einstein conden- 
sate accompanied by its thermal cloud. Jackson and Adams have proposed to combine 
the GPE with a quantum version of the DSMC to numerically evaluate the dynamics of 
such a fluid. Numerical studies based on a generalized GPE combined with a semiclassical 
kinetic equation, and including collisions between the condensate and the thermal cloud, are 
also becoming avaible for dynamical properties of a trapped Bose-condensed gas at finite 



temperature 



In the present work we proceed along the path traced by Jackson and Adams [^. We 
propose a different approach to the solution of the GPE and a different method for preparing 
the initial equilibrium state, which would be immediately applicable to a multi-component 
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cold-atom cloud. The method is applied to two classes of dynamical problems: the first 
concerns the ballistic expansion of a fermion cloud and the role played by the presence of 
a condensate, while the second concerns a Bose-condensed gas in a periodic optical-lattice 
potential at finite temperature. 

After introducing in Sec. |I| the model for both the equilibrium state and the dynamical 
evolution of the fluid, we describe in Sec. ^ the numerical methods that we have used to 
consistently solve the GPE for the condensate and the VLE for the Wigner distribution 
function of the cloud. The procedure followed in the actual computations is also outlined in 
Sec. |T|, with due emphasis on the preparation of the initial equilibrium state. The physical 
applications that we have carried out are presented and discussed in Sec. A discussion 
of computational aspects in Sec. M and some final remarks in Sec. conclude the paper. 



II. THE MODEL: A REVIEW OF THE MAIN CONCEPTS 



A. Transport in a normal quantum fluid 



In the original formulation by Kadanoff and Baym (KB), the problem of transport in 
a quantum fluids in the normal state is tackled by deriving an analogue of the Boltzmann 
equation for the Wigner distribution function /(p, r, t) from the microscopic equations of 
motion for the non-equilibrium density matrix p(ri, ti; ri'ti') =< ^/'^(ri, ti)'?/'^(ri/, ti') >, 
which is defined through the particle creation and destruction operators ip^^ in the presence 
of an external, slowly varying disturbance U{r,t). Namely, 

/(p,r,t) = y"dxexp(-2p-x) <^i(r + x/2,0)^^,(r-x/2,t) > (1) 

where r = (ri + ri/)/2 and x = ri — ri' are the center-of-mass and the relative coordi- 
nate of the two particles. The moments of the Wigner function yield observables such 
as the particle density n{r,t) = (27r)~^ J dpf{p,r,t) and the current density j(r,t) = 
(27r)-3/dp(p/m)/(p,r,t). 

Contact with the Boltzmann transport equation is made by performing gradient ex- 
pansions. As in the conventional Boltzmann-equation approach, the validity of the KB 
formulation is limited to slowly varying perturbations. On the other hand, the advantage 
of the KB formulation is that higher-order correlations enter the equation of motion for the 
density matrix in an explicit manner, and therefore systematically improved approximations 
which are consistent with the conservation laws are accessible. Examples of such treatments 
of the correlation term are the Hartree approximation and the Born-collision approximation. 
In the former case the collisionless Boltzmann equation is recovered, while in the latter the 
collisional Boltzmann equation is extended to non-dilute systems by including the effect of 
the external potential on the motion of the particles between collisions. 
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B. Extension to coupled condensate-noncondensate dynamics 



The extension of the KB treatment to gases including a Bose-condensed component has 



been made by Kane and Kadanoff |T9| and further developed by Griffin and coworkers p0|j21 



and also by Wachter et al. p2[ through a different derivation. The presence of two com- 
ponents and the appearence of off-diagonal elements in the density matrix (the so-called 
anomalous densities) below the Bose-Einstein condensation temperature requires the intro- 
duction of three Wigner distribution functions: /c(p, r, t) for the condensate component 
described by < i})^'^ > and involving \ < ip^ > /f,(p,r,t) for the noncondensate described 



by the fluctuations operators ip^^ = tp^'— < > and involving < ipljtpu >, and /m(p, r, t) 
for the anomalous part involving < ipu'^u > and its Hermitean conjugate. 

We thus have to deal with the density of condensate nc(r, t) = {2n)^^ J rfp /c(p, r, t), 
the density of noncondensate nb{r,t) = {2n)^^ J dp fb{p,r,t), and the anomalous density 
m{r,t) = {27i)^^ f dp fm{p,r,t). Analogous expressions holds for the current densities. 

As to the consistency of the approximations with the conservation laws, the same general 
remarks as for normal systems apply. However, the appearence of the condensate introduces 
an additional principle of gauge invariance, leading to the requirement that the excitation 
spectrum to be gapless [^. It is well known p0|j25[| that approximations capable of simul- 
taneously accommodating the conservation laws and the gaplessness condition are hardly 
available, so that a choice has to be made depending on the specific conditions of density 
and temperature of the system. 

In the regime that we address in the present work the anomalous densities can be ne- 



',(t) 



',(t) 



glected, resulting in the gapless Hartree-Fock-Bogoliubov-Popov approximation [^. Thus, 
the equation of motion for the condensate wavefunction $(r,t) =< ipyir^t) > is 



ih 



9$(r,t) 
dt 



2m 



v,t) 



<l>(r,t) 



(2) 



and is coupled to the coUisionless Vlasov equation for the noncondensate Wigner function 
fbip,r,t), 



+ ^ ■ Vr/,(p, r, t) - V^V^f^ir, t) ■ Vp/,(p, r, t) = . 
ot m 

The mean-field potentials in Eqs. (H) and (Rp are 



Vr\r) + W(r, t) + t/,K(r, t) + 2n,(r, t)] 



(3) 



(4) 



and 



V,'"{r, t) = Vr\r) + W(r, t) + 2Ug K(r, t) + n,(r, t)] 



(5) 
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including the time-dependent driving potential L{{T,t) and an axially symmetric confining 
potential given for a harmonic trap by V^^^{r) = nibuKr'^ + elz^)/2. In Eqs. @ and @, 
nc(r, t) = |$(r, t)p and Ug = 47r/i^a{,{,/mf, is the boson-boson interaction parameter in terms 
of the s-wave scattering length abb and of the boson mass m^. 

Once the algorithm to solve Eqs. (H) and @ is implemented, it is easily extended to a 
mixture of condensed bosons and a fermionic cloud in the coUisionless regime. In this case, 
the effective mean-field potentials become 

V:yir, t) = Vr\r) + W(r, t) + f/,n,(r, t) + Ufnf{r, t) (6) 

and 

V;"ir,t) = Vf{v)+Uf{r,t) + Ufn,{v,t) . (7) 

Here, V^^\r) = mfUjj{rj_ + ejz'^)/2 and W/(r, t) are the external trapping and driving poten- 
tials acting on the fermions and Uj = 27Th'^abf/mr is the boson-fermion coupling constant 
with abf the boson-fermion s-wave scattering length and rrir = mbmf/{mb + ^/), being 
the fermion atomic mass. Notice that fermion-fermion collisions in the s-wave channel are 
effectively suppressed by the Pauli principle in a dilute gas of spin-polarized fermions, as is 
relevant to current experiments on boson-fermion mixtures. 



C. Validity of the model 

To summarize, Eqs. (H) and (|]) describe the coupled dynamics of a Bose-Einstein con- 
densate and a bosonic or fermionic cold-atom cloud. In the former case the potentials V^l^ 
and V^^^ are used, and in the latter these are replaced by V^^^^ and V^^. The range of valid- 
ity of this approach is in principle limited to: (i) slowly varying space- and time-dependent 
external potentials, which allows a low-order gradient expansion of the equations of motion 
for the one-body density matrix; (ii) a coUisionless regime, which allows expansion of the 
self-energies to first order in the strength of the atomic interactions; and (iii) not too low 
temperature, so that the anomalous averages can be neglected. 

We implement below a numerical method to solve Eqs. (0) and (|^) and test it against 
dynamical behaviors in a one-component fermionic cloud and in clouds of either fermions 
or thermal bosons accompanied by a condensate. Of course, the role of the statistics enters 
at this level only from the initial distribution of the particles in phase space. We thus turn 
to discuss the equations that we use to determine the initial conditions for the subsequent 
time evolution. 
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D. Semiclassical description of the equilibrium state 



Before proceeding to present the algorithm for our dynamical simulations we briefly 
recall the basic steps that we take in preparing the initial state of the gas in thermodynamic 
equilibrium and in evaluating the corresponding densities of the condensate and of the 
fermionic or bosonic cloud. We refer the reader to Ref. p6| , ^ for the details of the theory 
and for a discussion of the excellent agreement that it yelds with thermodynamic data on 
Bose-Einstein condensed gases, under conditions of temperature and dilution that will be 
verified in the calculations reported in Section IV below. 

The equilibrium condensate density is calculated within the Thomas-Fermi approxima- 
tion, which amounts to neglecting the kinetic energy term in the GPE. Its validity is ensured 
whenever the average mean-field energy UgUc is much larger than the typical confining en- 
ergy. In the case of harmonic confinement the condition Ncatb/dho ^ 1 is required, where 
Nc is the number of atoms in the condensate, a^o = {h/mbtObY^'^ is the harmonic oscillator 
length and ZJb = ^^b^^^^ is the geometric average of the trap frequencies. The equilibrium 
density profile of the condensate is given by 

n,(r) = ^[f^b- Vr\r) - hjUbjir)] ^(/i, - Vr\r) - hjUbjir)) (8) 

where kb = '^Ug, kf = Uf and fib is the chemical potential for the bosons. 

For the equilibrium cloud density of bosons or fermions we adopt the semiclassical 
Hartree-Fock scheme [^]. This choice is justified as long as the gas is in a very dilute 
regime, a condition which is usually met in current experiments. In this approximation we 
have 



Tl} , (9) 



with V^'^'^-^(r) and Vj-^-^{r) determined by the confining potentials supplemented by static 
mean-field interaction term as in Eqs. (|5|) or (|^). 

The chemical potentials fib and fif are determined from the total numbers of bosons and 
fermions. In the case of a bosonic thermal cloud, fib is fixed by the relation 

Nb = j drK(r)+n,(r)] . (10) 

For a fermionic cloud the chemical potential is determined from the total number of fermions, 

Nj = J dvNfir). (11) 

These equations complete the self-consistent closure of the model in the initial equilibrium 
state. 
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III. THE NUMERICAL METHOD 



In this section we present the numerical procedure that we have used to solve the sys- 
tem of Vlasov-Landau and Gross-Pitaevskii equations. Since most experimental setups are 
invariant under rotation in the azimuthal plane, we use cylindrical coordinates {r,z}. The 
wavefunction $ is discretized on a two-dimensional grid of A^^. x A^^ points, which are uni- 
formly distributed in a box of sizes r^ax x '^Zmax- that is, = $(rj, Zk) with 

rj = ij-l)Ar {j = l,...,Nr) ^^^^ 

Zk = -Zmax + {k- l)Az (k = 1, . . . , N^) 

Ar and Az being the steps in the two space variables. The particle distributions are dis- 
cretized by means of a set of P computational particles, 

p 

/(p, r, t) ^ /p(p, r,t)^Yl - W)5(P " P. W) (13) 

1=1 

where the 6P phase-space coordinates rj(t) and Pi(t) represent the actual numerical un- 
knowns entering the VLE. They obey Newton's equations, 

— = Pi m 

dt ' 

where Fj collects all forces acting on the i-th computational particle. As a result, at each 
time step we have to solve for Nq = N^Nz grid unknowns coupled to 6P discrete particle 
coordinates. In order to keep the systematic signal reasonably above the noise level, each grid 
cell contains of the order of 10—100 computational particles. As a result, P ~ IOA'g — IOOA'g- 



A. Propagation of the Gross-Pitaevskii equation 

The GPE is advanced in time by an explicit finite-difference method based on a non- 
staggered variant of the Visscher method The full details are given in original papers 
1^,|T0| and here we shall just review the main ingredients of the algorithm. The basic idea 



is to advance the real and imaginary parts of the wavefunction, say A and B, in alternating 
steps as follows: 

BJ,^'-BJ-' = 2[K,,iA) + VJl{A)]At ^ ^ 

where n = 1, . . . , labels the discrete time sequence t„ = nAt. The scheme is initiated as 
follows. The initial conditions specify the values of A and B at n = and subsequently a 
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first-order Euler step provides their values at n = 1. With these values available, all later 
steps n = 2, 3, . . . , Nt are taken by using the above set of equations. 
In Eq. (|T3|) Kjk is the kinetic energy operator, 

K (^\- r '^J-i.fc I '^i+i'fe I "^i-fe-i I "^i-fc+i or ^ I M * 1 (^(\\ 
Zm Ar Ar Az Az Ar Az 

and 

V,km = V'ff{r^,Zk;njk)<^,k (17) 

is the potential energy operator, including both external and self-consistent interaction 
terms. The self-consistent potential requires the specification of the particle density njk 
at each node of the spatial grid. This is obtained by convoluting the discrete particle distri- 
bution, 



njk 



W = E W^i'Mt) (18) 



where Cjk is the grid cell centered in rj+1/2 = {rj + rj+i)/2, and ,2^+1/2 = + Zk+i)/2, 
while /j = 1 if the particle belongs to Cjk and fi = otherwise. The factor Wjk,i weights 
the contribution of particle i to the density at the grid-point {rj^Zk). 

We adopt a bilinear " cloud-in-cell" (CIC) interpolator (see fig. [l|a), which yields the 
scattering (particle-to-grid) rule p3 



Wjk,i = ei{rj)ei{zk) {j=ji,ji + l, k = ki,ki + l). (19) 

Here (jj, fcj) identifies the lower- left corner of the grid cell to which the i-th particle belongs. 
In Eq. ( p!UD ej(r) and ei{z) are one-dimensional piecewise linear splines centered on the 
discrete particle position (rj,2j), 

1 + — — {ri - Ar < r < r^) 



Ar 

A 

ri < r < ri + Ar) 



Ar 



and similarly for ei{z). 



B. Propagating the Vlasov-Landau equation 



The VLE is advanced in time by a standard Particle-in-Cell (PIC) method ||2^. In 
a modified Verlet time-marching scheme, we obtain the following set of discrete algebraic 
equations: 
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(21) 



' ri{t + At) = riif) + v„{t)At + ari{t)Af/2 
{riei){t + At) = (ri^i)(t) + veiit)At + aei{t)At'/2 
Zi{t + At) = Zi{t) + v,i{t)At + a,i{t)At'/2 
Vriit + At) = v„{t) + a„{t)At 
veiit + At) = ve^it) + aei{t)At 
v,i{t + At) = v,i{t) + a,i{t)At 

Here ar,a0, az and f ,., fe, Vz are the three components of the acceleration and of the velocity 
along the (r, 9, z) coordinates. The advantage of the modified Verlet time marching is 
that fourth-order accuracy is preserved while synchronously keeping the coordinate and 
momentum degrees of freedom on the same sequence of discrete times. 

The algorithm is standard except for the specification of the self-consistent coupling, 
namely the force due to the density gradients. In the first place, we form density gradients 
from auxiliary values of the density field %+i/2,A:+i/2 at cell centers. 



Gz 



[^i+l/2,fc-l/2 — "'j-l/2,fc-l/2 + %+l/2,fc+l/2 — "-j-l/2,fc+l/2]/ (2Ar) 
["-i-l/2,/c+l/2 — ^^^-1/2,^-1/2 + %+l/2,A:+l/2 " ^j-l/2,fc-l/2] /(2 A^:) 



(22) 



The azimuthal acceleration is zero in cylindrical symmmetry. Next, the grid-forces are 
evaluated on the discrete particle locations, this being the inverse of the scattering operation 
discussed in the previous section (see fig. [^). The grid-to-particle convolution is 



i,{j+s,k+s)^ j+s,k+s- 



(23) 



In order to avoid spurious self-forces we use again CIC-interpolation, which amounts to using 
the same weighting function as for the GPE in the grid-to-particle scattering rule, 

Wijk = ej{ri)ek{zi). (24) 

With the force/acceleration field transferred to the particle locations, everything is set to 
march the VLE in time. 



C. Boundary conditions 

The conditions imposed on the wavefunction are (i) periodicity along the z coordinate, 
(ii) symmetry (zero radial gradient) at r = 0, and (iii) vanishing at the outer radial boundary 
r = R. For the discrete particle distribution we have again periodicity along z and specular 
reflection at the outer radial boundary. Specular reflection means that a particle flying from, 
say, Ta < R to r/3 > R is replaced by a particle at r = 2R — r/3 with inverted radial speed. 
Since the particle trajectories are tracked in a three-dimensional cylindrical coordinate frame, 
the r = axis requires no special treatment. 
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1. Time-step considerations 



The GPE and the VLE are advanced on the same discrete time sequence. This maximizes 
simphcity, but imphes that the time-step is controlled by the fastest process at work, which 
usually is the self-consistent potential acting upon the condensate wavefunction. Better 
efficiency can be achieved by sub- cycling the time-stepper, namely by advancing the slowest 
equation (say the VLE) only every Aty/^/At^p steps, Atyi and Atcp being the largest time- 
steps allowed by the stability conditions on the two equations. The maximum time-step for 
the GPE solver is estimated from 

AtGp{C^^ + ChX^) < 1 (25) 

where 5 is a typical mesh size, Vm is the maximum value of the potential and Ci and C2 
are two coefficients 0(1) which depend on geometry and dimensionality. We note here the 
concurrent effects of quantum diffusion (kinetic energy) and scattering/absorption (potential 
energy) . 

The maximum time-step for the VLE solver is estimated from 

AWl'^ < 1 (26) 
d 

where Vmax is the maximum speed in the velocity grid, of the order of the Fermi velocity 
for fermions and of the thermal velocity for bosons. Under ordinary conditions the kinetic 
energy contributions dominate over those from the potential energy, so that the condition 
for the GPE to be the time-limiting section of the code takes the form of a numerical 
"uncertainty principle" , mVmax^ > ^- For the cases discussed in this work, this inequality is 
generally fulfilled within a factor ten, so that sub-cycling is not compulsory. The inclusion 
of coUisional interactions would make it mandatory. 



2. Radial singularity 

A source of potential trouble is the singularity in r = 0, which is known to affect all 
calculations in polar coordinates. To date, singular factors 1/r arc regularized by a simple 
numerical cutoff 1/r ^ l/(r + rc) with ~ O.OOlAr, with a check that the physical results 
are virtually insensitive to the specific value of Vc- 

Another undesirable sidc-cffcct of the cylindrical geometry is the relative depletion of 
near-axis cells, which tend to host fewer particles just because of the rAr volumetric effect. 
On the other hand, this volumetric depletion is often more than compensated by the physical 
behavior of the radial density, which is generally largest at r = 0. At any rate, the volumetric 
effect can be readily disposed of by moving to a non-uniform grid along the radial coordinate. 
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3. Statistical noise 

Considerations of statistical accuracy require of the order of a few tens of particles per 
grid-point (or equivalently per grid-cell) in order to keep the noise-to-signal ratio below an 
acceptable threshold. A practical consequence of this statistical accuracy requirement is that 
the VLE part of the computational scheme should be designed in such a way as to evolve 
these tens of computational particles in approximately the same amount of CPU time that 
it takes the GPE solver to advance a single grid-point. 

Another interesting consequence is that - at variance from ordinary situations in (clas- 
sical) rarefied-gas dynamics - the number of computational particles in the simulation of 
Bose-Einstein condensates far exceeds the number of physical atoms, typically by a factor 
10^ in our case. As a result, each single computer simulation performs de facto a built-in 
ensemble average over a set of about a thousand realizations. 



D. Procedure: preparing the initial state 

The initial condition for the populations of bosons or fermions in the cold-atom cloud is 
prepared as follows. Particles are sampled from the probability distribution functions (pdf ) 

where i] = Pifibj — Veff +j9^/2m) with (3 = l/fc^T the inverse temperature. 

The initializiation procedure starts by assigning to each spatial cell centered about po- 
sition r a corresponding amount of computational particles, 

AiV(r) = nbj{r)AV{r) (28) 

where AV{r) = 27irArAz is the cell volume. These particles must be sampled in momentum 
space according to the pdf (p7|). Owing to the non-separability of the pdf, straightforward 
sampling based on exact inversion is ruled out and more general - but less efficient - ac- 
cept/reject methods must be resorted to. 

The particle momenta are sampled from the distribution function using the standard 
Box- Mueller algorithm in three dimensions for cylindrical coordinates [Q, 

P± = Pmaxn 

e = 27rr2 (29) 

Pz = PmaxC^rs - 1) 

where ri, r2 and are sampled from a uniform distribution in the range [0, 1]. The maximum 
momentum Pmax is taken to be of order 2pF for cold fermions and of order SrribVth for bosons. 
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with pp = hkf being the Fermi momentum and Vth the average thermal velocity. The particle 
momentum coordinates in the azimuthal plane are evaluated as 

Pr = p±sin^ ^^^^ 
Pe = p±cos^. 

Finally, in order to avoid poor acceptance rates, the standard accept/reject test is per- 
formed by comparing the pdf with the maximum value /max(r) that it can take in the cell 
at position r, that is 

If f{Pr,Pe,Pz]r) > r4fmax{r) : accept 

Else : reject (31) 

where r4 is a random number uniformly distributed in [0, 1]. 



IV. PHYSICAL APPLICATIONS 
A. Expansion of cold fermions 

As a first test of the dynamical algorithm, we consider the expansion of a cloud of 
cold fermions after release of the harmonic trap. The neglect of the coUisional integrals is 
a reliable approximation in this context, since as already remarked the s-wave scattering 
between spin-polarized fermions is suppressed by the Pauli principle. 

This problem has been analytically solved in Ref. [^] under the assumption of ballistic 



expansion. The time-evolution of the mean square radii is found to be 

<rHt)>=±-E,.,^^(l+.]e) (32) 

and 

< z^t) >= ±.E,..-^^il + i,,.,)V) . (33) 

In Eqs. (^2]) and ( |33D Erei is the so-called release energy, namely the energy of the system 
with Nf fermions after switching off the trap, which amounts to one half of the total average 
energy. For a non-interacting fermion gas at temperature T > 0.2Tp this is best approxi- 
mated by the classical relation E^ei — "iN fksT /2, while at lower temperature it is given by 
Erei ^ {?>/A)NfEF[l + (27r/3)(T/T^)2], with Ep = keTp = {QNffl^njjjf being the Fermi 
energy. 

We prepare the initial thermodynamic state (^) for Nf = 1000 '^^K atoms in an isotropic 
trap with c<j/ = 27r x 15.92 rad/s at T = 0.55 Tp ^ 7.6 nK. The chemical potential of the 
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gas is /i/ = 0.063 Ep = 1.14 fkuf. We use &NrX mesh of 201x401 points with 1.6 ■ 10^ 
representative particles in a box of sizes 40 and 80 in units of aho = \Jh/ {nifUOf) along 
the radial and axial directions, respectively. The time-step in the dynamical simulation is 
ujf/^t = 10"^ 

We then evaluate from the simulation runs the radial width 



arAt) = \l j dxdydz[{x{t)~ <x{t)>f + {y{t)- <y{t)>f]nf{v-t) (34) 
and the axial width 



a,{t) = \j j dxdydz{z{t)- <z{t)>ynf{r;t) (35) 

of the cloud as functions of time, after averaging over the density distribution nf{r; t). Here, 
the center-of-mass coordinates are defined as < x{t) >= J dxdydzxnf{r]t) and similarly 
for <y{t) > and < z{t) >. Of course, during free expansion the center-of-mass coordinates 
must remain unchanged: this property is used as a test of the numerical method. 

Fig. I shows that the calculated (circles) and (squares) agree with the results 
from the analytical expressions ( P^ and ( |5BD (solid lines), where the classical expression has 
been used for Erei- Snapshots of the density profiles at selected times are shown as contour 
plots in fig. |. The definition of the profile degrades in time because the number of particles 
per cell drops during the expansion. 

After this test of the reliability of the simulational method, we proceed to use it in some 
original applications. 



B. Expansion of a mixture of a condensate and a cold-fermion cloud 

As a first novel application we look at the case in which a core of Bose-condensed atoms 
is present inside the dilute Fermi gas. We prepare a state with Nf = 1000 '^^K atoms 
and Nc = 10^ ^^K atoms in identical harmonic traps and at the same temperature as for 
the Fermi gas studied in Sec. IV. A. The scattering length which describes the interactions 
between the atoms in the condensate is abb = 80 aBohr, while the interspecies scattering 
length is ttbf = 40 asohr- 

The inclusion of the GPE algorithm at fixed mesh size normally requires shorter time- 
steps to mantain stability. We make the choice of a thinner 501x1001 mesh than in Sec. IV.A 
in order to keep the time-step at uAt = 10~^, all other simulation parameters remaining 
the same. 

The initial state is characterized by fi^ = 0.52Ef = 9.51 hwf and fij = O.lOEp = 
1.83 hjf. We display in fig. ||the average widths (circles) and o"^ (squares) for both the 
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fermionic cloud (open symbols) and the condensate (filled symbols). The solid lines are the 
analytical solution for the ideal Fermi cloud as in fig. It is seen that with the scattering 
lengths of the ^^K-^°K mixture the mean-field force of the inner condensate core on the outer 
fermionic cloud is not strong enough to sizeably affect the expansion of the latter. 

Snapshots of the condensate density profiles and of the fermionic cloud at times t = 0, 
8.5, 17 and 25.5 ms are displayed as contour plots in figs. ^ and ^. Comparison with those in 
fig. I shows that the reduced number of computational particles per cell tends to increase the 
statistical noise. This degradation worsens as the simulational time elapses, as is evidenced 
by the last snapshot in fig. ^ 



C. Oscillations of Bose gases inside an optical lattice 



Here and in the following subsection we apply our numerical method to study the dynam- 
ics of a Bose-Einstein condensate and a thermal cloud of ^''Rb atoms at finite temperature 
inside a one-dimensional optical lattice. The initial state is prepared by adding to the har- 
monic trap, described by l^rap(r) = (l/2)mco'^(r^ + e^^;^), a periodic potential given by 
Viatt{z) = tt-^ij sin^(/c£,2), where Eji = h^k\/{2m) is the recoil energy and = 27r/A is the 
wave number of the laser beam which creates an optical lattice with period Tr/ki in the axial 
direction. 

Such a system, which has been realized at LENS and examined numerically at T = 
by two of us shows a rich variety of dynamical behaviors. Thus, the study of the 

sloshing-mode oscillations of an almost pure condensate with N = 3 ■ 10^ atoms in a lattice 
with a = 1.6 shows that superfluidity is superseded by dissipation as the initial displacement 
of the condensate away from the harmonic-trap center is increased. This behavior is quan- 
titatively understood as a gradual destruction of superfluidity via emission of sound waves 
in the periodically modulated inhomogeneous medium Below the dissipative thresh- 

old, on the other hand, the oscillatory motion of the condensate through the optical lattice 
can be mapped into the dynamics of superconducting carriers through a weak-link Joseph- 
son junction [Q. This implies the possibility of observable resonances and of multimode 
behavior. 

Here we extend the above numerical studies to contrast the oscillations of a condensate 
with the motions of a thermal cloud. We prepare initial states for the two cases a = 1 and 
a = 5, both at T = for the BEC |Q and for the thermal cloud at temperature T above the 
critical temperature T^. We give an initial displacement Az = 42.6 fim to the trap center 
and follow the subsequent dynamics with a time-step uAt of order 10~^. 

The snapshots of the atomic density show that for a = 1 (see Figs. |^ and ^ the 
condensate behaves as a superfluid executing harmonic oscillations at a frequency equal to 
the trap frequency, while the thermal cloud at T > Tc diffuses away in a quarter of a period. 
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For a = 5 (see Figs. |^ and |TOD the condensate breaks instead into fragments as it attempts 
to perform the first oscillation, and after a period its center of mass becomes localized at 
the bottom of the harmonic well. In the same setup the thermal cloud becomes localized at 
the center of the trap in one tenth of a period and spreads out. 

Figure ^ gives a clear picture of these behaviors by reporting the axial center position 
and width of the condensate and of the thermal cloud as functions of time in the two cases. 



D. Expansion of a Bose-condensed gas in an optical lattice 

In our last study we look at the expansion of a Bose-Einstein condensate and its thermal 
cloud, which initially reside in a harmonic well and a superposed optical-lattice poten- 
tial. The external potentials are characterized by parameters typical of an experiment at 
LENS H, namely = 27r ■ 90 rad/s, e = 8.9/90, 27?/^^ = 795 nm and a = 5. The 



condensate contains 6935 ^^Rb atoms and the thermal cloud is composed of 3065 atoms: 
the temperature of the gas is T = 86 nK = 0.24 E^/kB and its chemical potential is 
fi = 5.86 fihj = 0.14 Er. We use a mesh of 111x2801 points with 308000 representative 
particles. 

We evolve the gas with a time-step uAt = 7 ■ 10^^ after switching off both the harmonic 
trap and the periodic potential. Snapshots of the atomic densities of the condensate and of 
its thermal cloud, taken at the moment in which the potentials are switched off and after 
3.5, 7 and 10.5 ms of free expansion, are shown in figs. O and |13l The condensate is seen 



in fig. |12] to develop side bands which separate out of the central cloud, while the thermal 
cloud in Fig. |l^ simply spreads out. These features of our numerical results reproduce those 
observed in the experiments [^ . 



The appearence of side bands in the condensate during expansion is due to Bragg scatter- 
ing against the periodic potential. In fact, in a long-time simulation run of a one-dimensional 
model of the expansion we have found that the condensate side-bands move at velocity 
V ~ 2/i/ci/m, corresponding to the momentum associated with the first reciprocal vector of 
the optical lattice. 



V. COMPUTATIONAL REMARKS 

We have assessed the computational performance of the numerical method by repeating 
the test of Sec. IV.A after changing either the number of computational particles or the 
mesh size. We list in Table 1 the computational times elapsed while running the HPF-PGI- 
compiled code on a fully dedicated IGHz Pentium III SCSI. 
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These data provide the following values for the specific CPU time costs of the GPE and 
VLE of the code per time-step: 



These figures invite a number of comments. First, they show that the VLE section can evolve 
just a few computational particles while a single grid-point of the GPE solver is advanced. 
Since statistical accuracy requires of the order of 10 — 50 particles per cell, we conclude that 
the VLE solver is a potential computational bottleneck. 

Let us nonetheless assume that the VLE and GPE sections can evolve on a one-to-one 
CPU time basis. We can then focus on the grid part only and estimate the feasibility of 
large-scale applications to finite-temperature condensates in optical lattices. Covering a 
simulation span of 100 ms in steps of 0.1 fis requires 10 timc-stcps. At a cost of 1 fis per 
time-step and grid-point, a grid with, say, 10^ grid-points takes of the order of 10^ seconds, 
namely almost two weeks of CPU time to complete. 

Ways to achieve substantial speed-up arc clearly needed. Among others, two promising 
(and not mutually exclusive) strategies are non-uniform meshes and parallel computing. 
Both strategies appear conceptually straightforward and will be the object of future work. 



The increasing complexity and variety of phenomena observed in current studies of the 
dynamical behavior of normal and superfiuid quantum gases at finite temperature motivate 
the development of suitable numerical tools to assist theoretical understanding. 

To this aim, we have combined a particle-in-cell method with an explicit time-marching 
algorithm to evaluate the time evolution in models of a Bose-Einstein condensate and a 
cold-atom cloud. 

We have tested the method against known analytical results in the simple physical situ- 
ations offered by the expansion of a coUisionless fermionic cloud without and with an inner 
Bose-condensed core. We have also applied it to simulate novel experimental observations 
on the dynamical behavior of a condensate with its thermal cloud in a harmonic plus optical 
lattice potential, where we have found substantial accord with current experiments. 

We have also analyzed those computational aspects of the algorithm which are most 
relevant to applications in large-scale problems. This analysis emphasizes the need for non- 
uniform meshes and parallel computing. On the physics front, an extension of the method 
to include the quantum coUisional integrals is under way. 




fis I grid — "point 
lis / particle 



(36) 



VI. CONCLUDING REMARKS 
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FIG. 1. Bilinear CIC interpolation to weight the contribution of particle i to the density at 
;rid-point {rj,Zk) (a) and to weight the effect of the forces at grid-point {rj,Zk) on particle i (b). 
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FIG. 2. Expansion of a cold fermionic cloud after release from a harmonic trap: radial and axial 
widths of the density distribution (in units of aho) as functions of time. Circles: 0"^^^; squares: (T^; 
solid lines: analytical expressions 
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FIG. 3. Expansion of a cold fermionic cloud after release from a harmonic trap: snapshots of 
the density distribution, shown as contour plots. Prom top left to bottom right: t = 0, 5, 10 and 
15 ms. The axial and radial coordinates are in units of a/jo- 
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FIG. 4. Expansion of a cold fermionic cloud and an inner condensate core after release from a 
harmonic trap: radial and axial widths of the density distributions (in units of o/io) as functions of 
time. Circles: Or^\ squares: cr^; open symbols: fermionic cloud; filled symbols: condensate; solid 
lines: analytical expressions (|3^ ) and (|3^ ) from Ref. |31|. 
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FIG. 5. Expansion of a cold fermionic cloud with an inner condensate core after release from a 
harmonic trap: snapshots of the condensate density distribution, shown as contour plots. Prom 
top left to bottom right: t = 0, 8.5, 17 and 25.5 ms. The axial and radial coordinates are in units 
of aho- 
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FIG. 6. Expansion of a cold fermionic cloud with an inner condensate core after release from a 
harmonic trap: snapshots of the fermionic density distribution, shown as contour plots. Prom top 
left to bottom right: t = 0, 8.5, 17 and 25.5 ms. The axial and radial coordinates are in units of 
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FIG. 7. Oscillations of a condensate in a harmonic plus shallow optical-lattice potential: snap- 
shots of the density distribution, shown as contour plots Prom top left to bottom right: t = 23.3, 
46.6, 69.9 and 93.2 ms. The axial and radial coordinates are in units of a^o- 
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FIG. 8. Oscillations of a bosonic thermal cloud in a harmonic plus shallow optical-lattice poten- 
tial: snapshots of the density distributions, shown as contour plots. From top left to bottom right: 
t = 6, 12, 18 and 24 ms. The axial and radial coordinates are in units of a/jg. 
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FIG. 9. Oscillations of a condensate in a harmonic plus deep optical-lattice potential: snapshots 
of the density distributions, shown as contour plots. Prom top left to bottom right: t = 23.3, 46.6, 
69.9 and 93.2 ms. The axial and radial coordinates are in units of aho- 
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FIG. 10. Oscillations of a bosonic thermal cloud in a harmonic plus deep optical-lattice potential: 
snapshots of the density distributions, shown as contour plots. From top left to bottom right: 
t = 2.8, 5.6, 8.4 and 11.2 ms. The axial and radial coordinates are in units of aho- 
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FIG. 11. Oscillations of a condensate and a bosonic thermal cloud in a harmonic trap plus 
optical-lattice potential with: axial center-of-mass coordinate and average axial width of the density 
distributions (in units of a/^o) as functions of time. Continuous line, condensate with a = 1; dashed 
line, condensate with a = 5; crosses for the thermal cloud with a = 1 and circles for the thermal 
cloud with a = 5. 
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FIG. 12. Expansion of a condensate interacting with its thermal cloud after release from 
harmonic plus deep optical-lattice potential: snapshots of the density distribution, shown as contour 
plots. From top left to bottom right: t = 0, 3.5, 7 and 11.5 ms. The axial and radial coordinates 
are in units of a^o- 
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FIG. 13. Expansion of a the bosonic thermal cloud interacting with a condensate after release 
from harmonic plus deep optical-lattice potentail: snapshots of the density distribution, shown as 
contour plots. Prom top left to bottom right: t = 0, 3.5, 7 and 11.5 ms. The axial and radial 
coordinates are in units of aho- 
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CPU-time (hh:mm:ss) 


1.6 • 10^ 


201x401 


8:35:27 


8- 10^ 


201x401 


5:33:55 


1.6 • 10^ 


401x801 


10:22:24 



TABLE I. Expansion of a cloud of fermionic atoms after release from the harmonic trap. 
CPU-time (third column) elapsed on a IGHz Pentium III SCSI for simulating 20000 time-steps, 
corresponding to 20 ms, for various numbers of computational particles (first column) and mesh 
sizes (second column). 
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